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Single molecule experiments on double stranded B-DNA stretching have revealed one or two struc- 
tural transitions, when increasing the external force. They are characterized by a sudden increase 
of DNA contour length and a decrease of the bending rigidity. The nature and the critical forces of 
these transitions depend on DNA base sequence, loading rate, salt conditions and temperature. It 
has been proposed that the first transition, at forces of 60-80 pN, is a transition from B to S-DNA, 
viewed as a stretched duplex DNA, while the second one, at stronger forces, is a strand peeling 
resulting in single stranded DNAs (ssDNA), similar to thermal denaturation. But due to experi- 
mental conditions these two transitions can overlap, for instance for poly(dA-dT). In an attempt 
to propose a coherent picture compatible with this variety of experimental observations, we derive 
analytical formula using a coupled discrete worm like chain-Ising model. Our model takes into ac- 
count bending rigidity, discreteness of the chain, linear and non-linear (for ssDNA) bond stretching. 
In the limit of zero force, this model simplifies into a coupled model already developed by us for 
studying thermal DNA melting, establishing a connexion with previous fitting parameter values for 
denaturation profiles. Our results are summarized as follows: (i) ssDNA is fitted, using an analytical 
formula, over a nanoNewton range with only three free parameters, the contour length, the bending 
modulus and the monomer size; (ii) a surprisingly good fit on this force range is possible only by 
choosing a monomer size of 0.2 nm, almost 4 times smaller than the ssDNA nucleobase length; (iii) 
mesoscopic models are not able to fit B to ssDNA (or S to ss) transitions; (iv) an analytical formula 
for fitting B to S transitions is derived in the strong force approximation and for long DNAs, which 
is in excellent agreement with exact transfer matrix calculations; (v) this formula fits perfectly well 
poly(dG-dC) and A-DNA force-extension curves with consistent parameter values; (vi) a coherent 
picture, where S to ssDNA transitions are much more sensitive to base-pair sequence than the B to 
S one, emerges. This relatively simple model might allow one to further study quantitatively the 
influence of salt concentration and base-pairing interactions on DNA force-induced transitions. 

PACS numbers: 82.37.Rs, 87.15. La, 87.15.A, 82.39.Pj 



I. INTRODUCTION 



In the recent decades, many experimental developments have been devoted to the manipulation and analysis of 
single molecules such as nucleic acids and proteins, proving to be invaluable tools to understand their statistical and 
mechanical properties pQ. They include atomic force microscopy (AFM) [2], fluorescence microscopy [3j, tethered 
particle motion [3] , as well as magnetic and optical tweezers [SHE] ■ Understanding the mechanical properties of nucleic 
acids at the nanometric scale is crucial because they play a role in both their biological functions and their packaging 
in a variety of circumstances, such as interaction with other macromolecules (e.g. histones or ribosomes), DNA 
cyclization, DNA looping in some genetic regulatory processes, or nucleic acid packing in viruses. Furthermore, single- 
molecule measurements on nucleic acids give the possibility to investigate their interactions with partner proteins [TJ[5]. 

The present work primarily focuses on the response of nucleic acids to an external force applied without torsional 
constraint, e.g., by optical tweezers or AFM. Such experiments on double-stranded DNA (dsDNA) at room temper- 
ature show a sharp, few picoNewtons wide, cooperative overstretching "transition" |57j at a given "critical' force of 
around 60-80 pN, accompanied by a sudden 70 % increase of the contour length [5J |B]. 

By applying force up to 800 pN, Rief et al. 0E] measured a second transition at stronger forces, which is hysteretic 
and consistent with a peeling of one strand from its complementary strand. They show moreover that the critical 
force of this second transition depends dramatically on the DNA sequence, from 150 pN for A-phage DNA to 320 pN 
for poly(dG-dC). They did not measure any 2d transition for poly(dA-dT). 

The nature of the first transition at 60-80 pN remains controversial because it is not definitely established whether 
it is a transition from the native B-form to a new form of unstacked DNA remaining in a duplex form (the so-called "S" 
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form for Stretched) , or double-strand separation leading to two single stranded DNA (ssDNA) similar to the thermal 
denaturation. This controversy [10] is reviewed in detail, e.g., in the Refs. [TTMMj . Some authors even argue that the 
nature of the transition depends on the loading rate, the slowest rates enabling equilibration and thus denaturation 
under force [THIIH]. In a very recent work [TTHT^j . two different transitions, both occurring at 60-80 pN for A-phage 
DNA, have been revealed experimentally, one is a hysteretic strand peeling whereas the other is a nonhysteretic 
transition that leads to S-DNA. The selection between this two transition depends on the DNA sequence and the salt 
concentration. Such study may thus reconcile the previous ones, once a careful comparison of both the experimental 
salt concentrations and the sequence of the studied DNAs will be done. 

Previous arguments for a B to S-DNA transition was that the part of the force-extension curve beyond the transition 
does not fit with a Worm-like Chain (WLC) model with the bending modulus and the monomer size of ssDNA [5j.f6j.l20J . 
However, the simpler stretching experiment on single ssDNA is already not easily fitted by the WLC model. Several 
explanations have been put forward, such as electrostatic interactions |21j . discreteness of the chain [22j . and non-linear 
bond stretching [53] . Using a variational continuous WLC model with linear stretching, Storm and Nelson [2U [55] 
were able to fit ssDNA for forces smaller than 0.1 nN but with an extremely low value for the monomer size between 
0.1 and 0.25 nm. 

Generally, the challenge is to propose a consensual mechanism compatible with this variety of experimental obser- 
vations. For instance, no adequate analytical formula exists for fitting dsDNA force-extension curves including the 
transitions. Since the experimental force-extension curves are measured for forces that range from to 800 nN, the 
model should scan both the low force regime where bending and entropic effects are important and the strong force 
regime where bond extension is non- negligible. Existing theoretical works use thermodynamical approaches [10) . in- 
terpolations [26j[27] between two (semi-)flexible chains using the Marko-Siggia interpolation formula [28], generalized 
Poland-Scheraga models with external force 129, !3U] , two-state continuous WLC models treated variationally [2H [55] , 
or three-state models [TU [50] . 

In the present paper, we first focus on ssDNA force-extension curves and provide an accurate analytical formula 
for fitting the stretching curves up to 1 nN (Section [n]). This formula includes bending rigidity (which plays a central 
role in the low force regime), discreteness of the chain, and non-linear stretching already studied by Hiigel et al. [53] 
(which has been shown to be very important for forces stronger than a few hundreds of pN). We then show that 
describing a hypothetic B to ssDNA transition with models that do not consider explicitly the increase of the number 
of sub-nucleobase degrees of freedom scanned in the high force regime is out of reach. 

Next, we use, in Section |lTlj a mesoscopic approach to model the B to S transition. This is a coupled Ising/Heisenberg 
model which is an extension of the mesoscopic model that we had introduced in 2007 for the description of temperature- 
induced melting of dsDNA j3T[[35]. This model is first solved exactly in Section III A using pseudo- analytic transfer 
matrix calculations, following the same formalism as Rahi et al. [33] . We then propose a simple analytical formula in 
the strong force approximation (Section IIIB), which is in excellent agreement with the exact results. 

Finally, we compare our formulas (the fitting procedure are summarized in the appendix) to experimental force- 
extension curves in Section IV and show surprisingly good agreement with several data sets for poly(dG-dC) and 
A-DNAs. Since our model is on the same footing as the one describing DNA thermal denaturation, we propose an 
analytical formula for the "coexistence" line in the temperature-force diagram. Our final remarks are given in the 
Conclusion. 



II. SINGLE STRAND DNA STRETCHING 



Before considering the transitions observed in dsDNA stretching experiments, we focus on ssDNA which is a good 
candidate for modeling a scmiflcxiblc chain under external load from to 1 nN. At very large forces, the discrete 
nature of the polymer is probed, and it has been shown that force-extension curves are satisfactorily modeled by a 
freely jointed chain (FJC) model in the strong force regime [12, 22, 53JIM1I35]. We thus focus on the discrete worm-like 
chain model (WLC), where the chain made of N monomers of size a is described by the effective Hamiltonian 

N-l 

Hdwlc = ^2 K b (l - tj ■ tj+i) - afU ■ z (1) 

i=l 

where tj is the normalized orientation of monomer i, and nt> is the bending elastic modulus. The applied force is 
taken to be along z, f = fz. The partition function is given by 

N-l r 

Z = Y[ a 2 dU e -^ DWLC (2) 

i— 1 
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The factor a 2 is the entropic contribution of the integration factor in phase space J d(rj_|_i — ri)<5(|rj_|_i — r$| — a)(. . .) 
a 2 / dt(. . .). The extension of the polymer along the direction of the force f is given by 



[rw - ri]) 



d\nZ 



The mean-squared end-to-end distance at zero force in the limit N — > oo reads 



(R 



a N 



1 + u(k) 
1 — u(k) 



(3) 
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where u(x) = coth(x) — 1/x is the Langevin function and n = {3nb is the adimensional bending modulus. The discrete 
persistence length is thus £ p = — a/lnu(/t). Since we consider a semi-flexible chain, we assume k ^ 1. 

At small forces, such that the chain is only slightly deformed in the z direction, the chain can be viewed as a linear 
string of Pincus blobs [5B] of size £ = l/((3f) with L > £ > l v > a, where L = aN is the contour length of the 
polymer. Hence in this regime where flanf < 1, the linear response theory is valid and yields for the extension 
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Note that this relation is modified in good solvent according to the scaling law z ~ L(f / f) 2 / 3 . In the following we 
neglect the eventual polymer swelling for DNA in water due to van der Waals (35] or electrostatic interactions [3"7] . 
Moreover the electrostatic contribution to the persistence length which might be important for flexible chains such as 
ssDNA is taken implicitly into account by choosing fitting parameter |38l I39j . 

At very large forces, / 3> 4Kb/a, we see from Eq. ([T]) that the bending energy can be neglected in the Hamiltonian 
and the freely jointed chain model is valid. The force probes the discrete nature of the chain since the Pincus blob 
size is much smaller than the monomer size a, £ <C k^T / {Ak^o) . One thus finds the classical Langevin result |40) 



z ~ L u(a/3f) ~ L 1 - 



a/3/ 



for / > 
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For intermediate forces, (h^T) 2 / (aKf,) < f < 4Kb/a, the two terms in Eq. |T]) should be taken into account, although 
the bending energy is negligible in the z direction. Following Marko and Siggia |28j . we use the approximation of large 

forces, which states that is mostly along z (\t iz \ ^> \t ix \, \U y \) and thus \t iz \ = y/l — t 2 ix — t 2 y ~ 1— (tj x +t 2 y )/2. By 

noting that k(1— t,-tj_|_i) = |«;(tj — t^+i) 2 , the partition function can be rewritten as Z — Jl^i 1 / dtjdtj+iP(tj, tj+i) 

where P(tj,t,-|_x) = a2 ^ F K{t x ^ : t Xt i + i)K(t y ^, t y .i+i) is the transfer operator to be diagonalized. The eigenvalue 
equation is 
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a e 



dt' x K{t x ,t' x )<t>x(t' x ) / dt' y K(t y ,t' y )ct> y (t' y ) = A(l> x {t x )ct> y (t y ) 



(7) 



where we define the adimensional force F = a/3f and 



K(t, t') = exp 



+ M2 F 



2 (t-t>r-j(t 2 +t> 2 ) 



(8) 



This transfer kernel problem has first been solved by Fixman and Kovac [41 using a mode decomposition for a discrete 
worm like chain with extensible bonds (without applied force) , and then extended to the force case in Ref. [HU [35] . 
For sake of clarity and as an introduction for the more complex case developed in the next section, where the mode 
decomposition is not applicable, we re-derive the calculation in the following by using a different approach, i.e. solving 
Eq. ^ directly in real space. 

We define (f>i(t) and <j>2{t') such that 



dt' exp 



r(*-0 



/\2 



F 2,f2 



(9) 



If ^(t) = exp(-a 1 t 2 /2) and <j> 2 {t') = exp(-a 2 t' 2 /2), we find 
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If F\ — F2, then cf>i — cf>2 is an eigenfunction for a% — — a = yj (F/2) 2 + kF and A = A 2 . The partition function 
in the limit -/V — > 00 is thus 
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Z ~ e 



NF 



2na 2 



(11) 



and the extension defined in Eq. ([3]) is 



z = L\ 1 — 



1 



v^ 2 + 4 K ,F 



(12) 



Note that Eq. (12) encompasses the very strong 



Again this formula is only valid for large forces, i.e. for F ^> 1/k. 
force regime F 3> 4k where the result of Eq. ^ is recovered. 

An interpolation formula for the whole range of forces can be obtained following Marko and Siggia |28j by inverting 
Eq. |12|, subtracting the first two terms of the small z expansion of F{z/L), and then adding the zero force limit 
Eq. |5|. This yields the discrete Marko-Siggia interpolation [53] : 
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4k 2 - y/l + 4k 2 



(13) 



Experimentally, one observes that, at very large forces, / > fc^T/a, the DNA starts to stretch elastically and the 
FJC result, Eq. Q, must be corrected so that z > L. Different models try to take this into account by using a linear 
correction [28 , or using an extensible DWLC model [34]. However, it has been shown recently by Hiigel et al. that, 
beyond elastic stretching, non-linear terms are necessary to fit force-extension curves of peptides or polyvinyl-amine 
as well as ssDNA [55]. A consistent way is to modify Eq. (12) according to 



z = L[l+ U nl (f)} 1 



1 



VF 2 + AkF 



where 



C/ nl (/) = 1.172777 / - 3.731836 f + 4.118249 f (/ in units of 10 nN) 



(14) 



(15) 



is extracted from Ref. 23J [by inverting their Eq. (2)] and is the result of ab-initio quantum calculations for ssDNA. 
Corrections are on the order of w 1% for / ~ 100 pN, of ss 10% for ~ 1 nN. 



To reconcile the low force regime where Eq. (13) is correct and the large force regime, Eq. (14), where non-linear 



elasticity is non- negligible, we fit the experimental data for ssDNA using the interpolation Eq. ( 13 ) with z replaced 
by z(l + U n i(f j) in the rhs. It yields: 
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(16) 



This non-linear equation is easily plotted using a parametric plot if we replace U n \(f) by ^(/dms), where /dms is 
defined in Eq. (13). We have checked that this approximation is extremely good for forces in the nanoNewton range. 



We use Eq. ( 16 ) to fit experimental data with only three fitting parameters, namely the polymer contour length L, 



the adimensional bending modulus k and the monomer size a. 

The comparison with experimental data taken from [23] (kindly supplied by R.R. Netz) is shown in Fig. [lja) and 
(b). The fit using Eq. (16) is quantitatively good (the curve remains within the experimental error bars) and yields 
the following values for the fitting parameters: L — 3.40 /im, k — 1.5 and a = 0.20 nm. This fit is very constrained 
since it is done on a very large range of forces, form to 1200 pN. 

A first remark is that, although the persistence length is quite small t p = — a/\nu(n) ~ 0.24 nm ps a, the role of k 
is non-negligible in the low force regime and setting k — (FJC model) leads to a poorer fit (data not shown). 

As shown in Fig.jljb), the interpolation (black curve), Eq. |l6| starts to deviate slightly from Eq. |l~3"l ) (red curve) 
for / > 150 pN. Simply put, as already said by Hiigel et al. |23j. the non-linear stretching starts to play a significant 
role. Moreover, for / > 400 pN, we have \j\jF 2 + AnF ~ F^ 1 < 0.05 which becomes smaller than J7 nl (400pA r ) = 0.04 
in Eq. (14). In other words, for / > 400 pN, the entropy becomes negligible, the influence of k and a on the fit is 



small and the extension-curve is dominated by the non-linear stretching. 
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FIG. 1: (a) Extension vs. force for a ssDNA. Data (symbols) are taken from Hiigel et al. 
to a fit using the discrete worm like chain interpolation with the non-linear bond elasticity, Eq. ( 16 1 (the red one corresponds 
to discrete Marko Siggia interpolation Eq. ( |13[ l). The parameters values are: L = 3.40 /im, k = 1.5, a = 0.20 nm. (b) Zoom of 
(a) for small forces, together with the strong force limit Eq. ( 12 l (green). 



Although the L and k values were expected, the effective monomer length a = 0.20 nm is much smaller than the 
distance between two consecutive bases in DNA, a ss ss 0.7 nm. We define a corrective factor b = a ss /a = 0.285 
which is the ratio between the expected ss base size a ss and the fitted monomer size value a. Note that Storm and 
Nelson [25] found similar values (a = 0.17 nm) by using a Ritz variational approximation and fitting only on the 
[0, 200 pN] range. Refs. (221 [23] focus on the non-linear elasticity which is significant for forces larger than 400 pN. 
By using a non-linear Freely Rotating Chain model at large forces they found a monomer size multiplied by 2 for 
polyvinylamin [55] H3] and a smaller monomer size by a factor 0.5 to 0.8 for peptides [35]. In our case, the actual 
monomer length probed by a strong applied force is thus 1/6 = 3.5 times smaller than the inter-base distance for 
ssDNA. In other words, the number of degrees of freedom N for such strong forces increases by a factor 3.5. 

This increase of N at large forces raises an important question about a transition B-DNA to ssDNA in stretching 
force experiments. How to model this change using a mesoscopic model? This increase probably occurs abruptly 
during the transition and is likely to include chemical modifications unaccessible to classical mechanics. 



III. ANALYTICAL MODEL FOR B TO S-DNA TRANSITIONS 



While the previous section casts some doubt upon the adequacy of a mesoscopic model to model the B to ssDNA 
transition, we show in this section that it is possible to describe the B to S transition using such type of model. We 
use a Ising-Heisenberg coupled model, which has already been used by us for the theory of DNA denaturation [3TJ [35] . 
Other works have used such types of models [24, 25 . In this model, each base-pair i is described by: (1) its normalized 
orientation tj with fij = (B^ipi) the solid angle with respect to a fixed reference frame (x,y,z), (2) its internal state 
<Ji = ±1 (corresponding to B or S base-pair internal state), and (3) its length a(<7j). By generalizing Eq. ([!]), the 
effective Hamiltonian is 

N-l 

j3'H=22 K(ai, <Ti+i)(l - tj ■ t i+ i) - Pa(ai)fti ■ z + Hi(ai, a i+1 ) (17) 

i=l 

Compared to Eq. (JlJ, the bending modulus nfa, Oi+i) now depends on the internal state of base pairs i and i + 1 
(we also note in the following tj • t i+1 = cos7i i i + i) and the additional term 

Hl((Ti, 0"i+l) = —J(Ti(Ti + l — ^((Jj + CTj+l) (18) 

is the internal Ising free energy associated with base pair i and its interaction with base pair i + 1 (2/i is the energy 
necessary to break one base-pair and 2 J is the energy of a domain wall) [3"T1 132j . 
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The transfer operator P is then defined by 

(Cli,<Ji\P\Q.i + i,a i+ i) = a 2 (a) cxp [-Hi(ai,a i+ x) + nfa, a i+ i)(cosji, i+ i - 1) + /3a(ai)f cos 0*] (19) 



= a 2 (a)exp 



-HM^i+i) + K{(Tl ^ l+l) (U-t l+1 ) 2 +f3a(a t )fU 



(20) 



A. Exact diagonalization of P 

The idea of dealing with the WLC model under forces in the spherical harmonics basis goes back to Marko and 
Siggia jlSj, even though their use in different but related fields of physics goes back to the 70s [53]. In the present 
case, we diagonalize P by using the decomposition of a plane wave in spherical waves: 



= KCOS 7;,i + l 



I oo £ 

J^Y, 1 ^^ 2 yu^+i)YLm (2i) 



which implies that Yp m is an eigenvector of e K 



cos 7i,«+i 



^i e «ccB 7M+1 y im(n .j = & i(„)y /m (n i+1 ). (22) 

Note that the prefactor of Yg, m in the rhs. term is the spherical Bessel function usually denoted by ii{k) and that we 
also used the notation e~ K i(_{n) = e ^ G *( K ) in Ref. [52"] . 

If / = 0, P is block-diagonal in each (£, m) subspace with matrix elements (we switch to lighter notations): 

(£m, a\P\£'m', a') = a 2 (a) exp [-G £ (k(ct, a')) -Ui{a, a')] (23) 

They depend on I but not on m. When diagonalizing each 2x2 block, the eigenvalues are denoted by Xg y ± and the 
eigenvectors by \£, ±). The partition function is Z = (in) N J2± + -O^i m case °f periodic boundary conditions 
and by Z — (4ir) N J2e ±(V\®> ^ z ) 2 ^e'± m case °f f ree en d s ( wnere 1^0 is the adequate free-end vector) [32]. Of course, 
boundary conditions are irrelevant at large N [45], and we have checked that finite size effects arc negligible for the 
DNAs studied in the experiments considered in the following. 

If / 7^ 0, P is block-diagonal, but blocks are now infinite because different values of £ are coupled. A block thus 



corresponds to a given value of m. Eq. ( 22 ) has to be adapted to this case [33J E2HH] (Note that contrary to [33] , we 



do not need to explicitly treat the two single strands in the B to S transition). In Dirac notations, we have: 
^l^riH^j^ = ^(-lr^Vt^ + 1)(2*' + l)e- K ie{n) £ (2* 1+ 1) ( * [ ^ * m ) i tl W 

where we have used Wigner 3j-symbols (and m < £,£'). Such use of Wigner 3j-symbols already appeared in Blume et 
al. [43] - In this expression, k means k(<7, a') and a means a(a). In practice, to diagonalize each infinite block, a cutoff 
^max ~ 8 on £ must be chosen (large values of £ do not give significantly different results). Once P is diagonalized, Z 
can be computed from which the extension z is derived. 



B. Strong force approximation 



Similarly to SectionjllJ we use now the description of P in tangent vectors t,. An extension of the spinor eigenvector 
equation is now, in a symmetric form: 



dt' 



a ^+J K ++ (t,t') a+a_e- J K+_(t,t') 
K_+(t,t') a 2 _e-»+ J K__(t,t') 



$+(t') 
*-(t') 



= A 



$+(t) 
*-(t) 



(25) 



where A is the eigenvalue and $<j(t) = 4> X o-(t)4>yo-(t) are the unknown eigenf unctions. The transfer operator 
K CT(J /(t,t') = K a( ji(t x ,i! x )K aa i{ty^t' y ) with a, a 1 =± is the generalization of Eq. (|8l where 

K aa ,{t,t') = cxp 



^(,y)(t-o 2 -^-Jw 2 



(26) 
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Searching for an exact diagonalization of the transfer matrix is difficult because a Gaussian wave such as in Section [TT] 
is no more an exact eigenfunction. This is due to the fact that cross-terms are not symmetric in t and t': Kx t -i{t, t') 7^ 
Ki.-i(t' ,t) [but Ki i(t,t') — K-i i(tf ,t)]. For the same reason, the Ritz variational scheme [MJ [25] does not work 
for this model with a Gaussian variational eigenfunction. 

Nevertheless, eigenfunctions are Gaussian in 3 limits: 1) the homogeneous chain (this has been proved in Section pi)), 
2) the zero force limit / = 0, and 3) in the freely jointe d chai n limit, = 0. Indeed, by inserting in Eq. (25 1 Gaussian 
wave functions (j) a = C a exp(— a a t 2 /2) and using Eqs. ( 9|10 ), we find 



where 
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t 2 /2 = Ae -a_t 2 /2 



(27) 

(28) 
(29) 



By assuming a±± = a± = y/ (F±/2) 2 + n±±F±, and dividing Eqs. ( 27|28 ) respectively by exp(— a±t 2 /2), one finds 

„(a + -a + _)t 2 /2 _ x 



C 4 



6*4 



2tt 
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; (a_-a_ + )t 2 /2 _j_ 



27T 



F- 
2 



= A 



(30) 
(31) 



It is then straightforward to check that in case (2) (/ = 0), we have a± T = = a± and in case (3) (n aa ' = 0) 

a± T = Fa J 2 = a± and the Gaussians of the cross-terms are equal to 1. 

To proceed further, we make the approximation that the Gaussians in Eqs. ( 30|31 ) are negligible and thus assume 
that they are equal to 1 for all the parameter values. This allows us to write an effective Ising model since Eqs. ptplj ) 
do not depend on t anymore. Hence, coming back to Eq. (25), the transfer matrix reduces to a simpler 2x2 Ising 
transfer matrix: 



r 1 e Ho+Jo e -J(>-& 



-Jq+5 p—fio+Jo 



(32) 



with force- and temperature-dependent Ising parameters (we switch to the subscripts B for B-DNA instead of +, and 
S instead of - for S-DNA): 
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(33) 
(34) 
(35) 
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(37) 



In the limit N — > 00, the partition function is then given by Z — A N where the largest eigenvalue of the Ising 
matrix is 58 



A = e r ° +Jo [ cosh 



1 A*o 



sinh fj,Q 



-4J 



(38) 



The "magnetization" and the two-point correlation of this effective Ising model are (see Ref. 
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cosh jUq + \/ sinh 2 /zq + e _4J ° 



where Eq. (39) yields the fraction of base-pairs in the B (resp. S) state 

<PB,s(F) = 

as a function of the force. The extension computed according to Eq. ^ is thus 



a B N 



2a B 



<Pb + 7 1 



— -) 



7 



2ob k B s + F/2 + a B 2a s k bs + lF /2 + a s 



(39) 
(40) 

(41) 
(42) 



which shows that the last term is only relevant close to the transition where ((TiCTi+i) ^ 1. 

Eq. (42 1 is the second important result of the paper. As it is constructed, this formula is an interpolation between 



several limits. First, the result for an homogeneous chain in the strong force approximation, Eq. (12), is recovered by 
setting k b = Kg — k bs an d 7 = 1 in Eq. (42). Second, in the zero force limit F = 0, Eqs. ( 33|36 ) simplify to 



Mo = M + - In 



~f 2 K B 



Jq = J 



1 



In 



"■SB 

kb us 



T = - hx{K 2 SB K B Hs) 



ln(27T7) 



5 = 



(43) 



which are the renormalized Ising parameters already found in |48j . The FJC model corresponds to Kb = «s — k bs = 0, 
and Eq. ( 42 ) reduces to 



a B N 



(fB + 7Vs) - p 



(44) 



which is Eq. ^ slightly modified to take into account to the two accessible values for the base pair length <Zj. 

Finally, far from the transition, defined as (ctj) = or cquivalently Mo(/c) = for infinitely long DNAs, that is for 
forces such that (<7j<7j-|_i) = 1, we find B- (or S-) stretching behaviour, for /j,q > (respectively fj-o < 0): 



1 



a B N 



a B N 



\JAkbF + F 2 
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V4^ s7 F + 7 2 F 2 



/ « fc 
f » fc 



(45) 
(46) 



This last result is identical to Eq. ( 12 ) provided that the extension and the force are renormalized by the S-monomer 
length as = Jclb- Eqs. ( 45|46 1 are a generalization of the result of Cizeau and Viovy [26] where the continuous 
Marko-Siggia interpolation, valid for range 1/k <C F <C 4k, was used. 

In experiments, large forces on the order of several hundreds of picoNewtons are applied to B-DNA such that helix 
stretching occurs. This stretching is related to the torsional elasticity of the double helix, such as for a helical spring. 
It is incorporated linearly, following Eq. (14), by replacing / by /(l + f/Es) [231 [2H] m the prefactors independent 
of t and t' in the matrix elements of Eq. (25), where Eb (in pN) and the related adimensional E B = Po-bEb is the 



stretching modulus in the B state, taken as a fitting parameter. Eq. (42) becomes [59 



a B N 
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Eqs. ( 42|47 l allow us to fit with a high accuracy the various experimental DNA stretching curves taken from the 
literature and also compare perfectly well with the semi-analytical exact formula, Eq. ( 24 ) , as illustrated in the next 
Section. 



47) 



IV. COMPARISON WITH EXPERIMENTAL FORCE-EXTENSION CURVES 



We now compare our theoretical approach to various experimental data. Some of them, coming from optical tweezers 
experiments, are extracted from [55]. Rief et al. [5] also conducted several experiments using AFM on A-phage DNA, 
poly(dG-dC) and poly(dA-dT) DNAs, in order to explore the role of base-sequence on DNA stretching. 
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FIG. 2: (a) Extension vs. force for a poly(dG-dC). Data (blue symbols) are taken from Rief et al. 2 . Solid curves correspond to 
the discrete Discrete Worm like chain interpolation, Eq. ( 13p , for B-DNA (red), S-DNA (green) and with non-linear extensibility, 
Eq. (16 1, for ssDNA (pink). The black curve corresponds to Eq. (47 1, where linear stretching is included as shown by the blue 
curve for pure B-DNA. The red symbols correspond to the semi-analytical calculation Eq. (24 1 with ^ max — 8. Parameters 
values are: L B = 0.14 /im, k b = 147, 71 = 1.89, k S = k B s = 3.8, 72 = 1.145, b = 0.285, k ss = 5.54/(2 7 i7 2 ) = 1.28, 
J = 1.7, Eb = 12 00 p N. Inset: Fraction of base-pairs in the S state vs. force, Eq. (41 1 and Ising correlation function 
f 1) defined in Eq. (40 1 (dashed curve), (b) Same as (a) for A-phage DNA. Data (blue symbols) are data taken from [2J. 
kb = 147, 7 = 1.88, k s = kbs = 4, /x = 3.85, J = 2.05, E B = 1400 pN. 



M = 4.5, 

1 — (<Ti<T, 

Parameters values are 



A. Overstretching transition for poly(dG-dC) and A-DNA around 60—80 pN 



First, we focus on the B to S transition which occurs for / = f c ~ 60 — 80 pN with small differences related to 
the DNA sequence. We compared our analytical result Eq. (47) to experiments made on polyGC taken from [5] [see 
Fig. |2ja)] and A-phage DNA [figs, gb) and Fig. and Eq.^47| leads to very good fits of experimental data. The 
fitting procedure is detailed in the appendix. 

To begin with, we have checked in Fig. [2ja) that the semi-analytical calculation using the result of the Section 
IIIA (red symbols) and the strong force approximation of Section IIIB (black solid curve) are superimposed. Note 
that the results of Section IIIA are done without linear elasticity for B-DNA. This is the reason why there is a 
slight difference before the transition since we have plotted only Eq. ( |47| a nd not Eq. ( 42 1 for sake of clarity. It 
proves the validity of the strong force approximation used to derive Eq. (47) for the B to S transition. This is due 
to the fact that the transition occurs in the force range (60-80 pN) where / 3> /b = k B T '/ '(obH-b) — 0-1 pN and 
/ > fs = kBT/(a S Ks) ^ 2 pN. 

This is also confirmed by the plot of (fis(f) m the inset of Fig.[2ja), where both results are superimposed. Moreover, 
the superimposition of Eq. (47), valid for N — > 00, and the results of Section IIIA computed for N ~ 700 provide the 
undisputed evidence that the very small correction due to finite N lies within error bars. Note that the transition 
is abrupt as shown by the plot of 1 — (crj<7j_|_i) in Fig. [2]ja), and Eq. (47) gives the two correct limits far from the 
transition Eqs. ( 45|46 ). In the fitting procedure, the parameter kbs plays a similar role as the parameter J (see 
Eq. (34)). This is the reason why we chose kbs = K s- 

Furthermore, we notice that the fits of Figs. [2] and [3] yield similar values for 7, between 1.72 and 1.89, which are 
also comparable to those obtained by Storm and Nelson S24, 25 a for A-DNA (1.7-1.8). While the bending modulus of 
B-DNA is taken to be 147 k B T, the S-DNA one is much smaller, between 3.8 and 4 k B T. Finally, similarly to [24j. [25] 
we find Eb ~ 1000 pN. Contrary to Refs. [Ml we do not need to introduce an linear modulus for the S-form, 
which can be attributed to the fact that a continuous WLC model was used in [Ml [23, instead of a discrete one. 

By fitting the transition using Eq. (47), one finds 3.85 < fi < 4.5 in k B T units, which are reasonable values compared 
to that extracted from the Poland-Scheraga model [49] and fits of denaturation curves [50]. It is roughly twice the 
value found for for poly(dA)-poly(dT) in [351 [55], and is consistent with the fact that GC base-pairing energy is larger 
than AT one. The cooperativity parameter J is 1.7 < J < 2.05 which is also reasonable and a little smaller than the 
value of 3.6 found for poly(dA)-poly(dT) in [35] 0S]. The small variation of the fitting parameter values from sample 
to sample are probably due to the slight differences in DNA sequences and salt conditions. In Fig. |4^a) is plotted, 
in the (T, /) plane, the coexistence line defined by setting fio(T m ,F c ) = in Eq. ( |33| with the parameters values of 
Fig. [5Ja). It shows the same behaviour as in Refs. [Ml [33], with a re-entrance for (unreachable) high temperatures, 
and decreases linearly in the accessible temperature window [19] , 
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FIG. 3: Extension vs. force for a A-phage DNA. Data (symbols) are due to Cui and Bustamante for (a) and Ouzel et al. [5] 
for (b), and both taken from Storm and Nelson [25] . The solid curves correspond to the discrete Discrete Worm like chain 
interpolation, Eq. (|16[ ), for B-DNA without stretching (red), with stretching (blue) and ssDNA (green). The black curve 
correspond to Eq. \47[ , where linear stretching is included. Parameters values: Kb = 147, ks = k bs = 4, l/i?s = 0, and (a) 
7 = 1.795, n = 4.015, J = 2, E B = 1300 pN; (b) 7 = 1.715, fx = 3.85, J = 1.9, E B = 880 pN. 



To conclude this Section, the fact that Eq. (47) allows us to fit the transition observed experimentally for poly(dG- 
dC) and A-DNA indicates that the second state is indeed a S state and not a ss state. Indeed, a good fit of a transition 
to a ssDNA state around 60-80 pN would impose a much smaller value of the monomer size as discussed in detail in 
Section [TTJ and below. 



B. Second transition for poly(dG-dC) around 350 pN 

Rief et al. 2 observed a second transition when they stretched a poly(dG-dC) at larger forces, around ~ 350 pN, 
as shown in Fig.Qb). They argued that this transition corresponds to a S to ssDNA transition, where the final state 
corresponds to one single ssDNA strand which remains tethered, the second strand being unpeeled [2]. Indeed, this 
second transition, which is very smooth between, roughly 200 and 400 pN, shows an hysteresis which varies with 
the applied pulling speed of the AFM tip. They also observed this second transition for A-DNA at a smaller force, 
f c ~ 150 P N. 

Following these arguments, we try to model this second transition. First we use the result of Section [TTJ where the 
value of the actual bond size is a = 0.20 nm. We then fit the experimental data at strong forces in Fig. [2]ja). One 
finds consistently L ss = 7172-^s = 2.\QLb and k ss = 1.23 ~ 5.54/(27172), where 71 = as/as and 72 = a ss /as and 
the value 5.54 is taken from [5TJ[52]. These two values corresponds to the distance between two adjacent base-pairs 
along the helix (ss 0.7 nm) and to the accepted bending modulus value of a single ssDNA strand (persistence length 
£ p ps 1 nm [B]). Note that the persistence length of ssDNA can vary a lot with the salt concentration [55]. 



Second, to fit the transition, we use Eq. (47) where the B state becomes the S one and the S state is the ss one. 



However, as explained in Section [TTJ the change of degrees of freedom from B to ssDNA should prevent the success of 
the fit. But, since the ssDNA form is obtained for / > 400 pN, following Section[nJ the entropy is not dominant for this 
force range, the stretching being essentially due to bond deformations modeled by the non-linear stretching Eq. (15). 
Hence, in the absence of any model with different degrees of freedom in S and ss states, it is reasonable, as a first 
attempt, to keep the same number of degrees of freedom for this case [see in Fig. [4](b)] : a — a ss — 7172&B = 2.16as 
(or 6=1). 

Within this hypot hesi s, we are able to fit approximatively the experimental curv e by replacing the linear stretching 
term, F/Eb, in Eq. (47) by the non-linear stretching one for ssDNA given by Eq. (15). The values of the fitting Ising 
parameters for this transition are fj, = 3.8 and J = 0. It indicates that this transition is not cooperative at all. This 
is consistent with the commonly accepted picture of a destacked S-DNA: during the S to ssDNA transition, only the 
breaking of the hydrogen bonds between base-pairs occurs, the aromatic rings being already destacked in the S state. 
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FIG. 4: (a) Phase diagram in the temperature-force space corresponding to the first transition of Fig. [2^a). (b) Fit (black 
curve) of the second transition for poly(dG-dC) using Eq. ( |47[ ) where the linear stretching term is replaced by the non-linear 
stretching one f or s sDNA given by Eq. (15 1. The green solid curve is the same as in Fig. [2ja) and the red one is the large 
force limit Eq. (12 1. The parameters are the same as in Fig. [2^ a), Nas = Nasji = 0.265 /im, ks = 3.8, 72 = 1.145, 
k ss = ksss = 1.28; and the fitting parameters are fj, — 3.8 and J = 0. Note that, compared to Fig. [2|a), the fit is poorer for 
ssDNA since no parameter b is introduced. 



C. Nature of the transition for poly(dA-dT) 

In Fig. [5] is displayed an attempt to fit the transition observed for poly(dA-dT). Following Rief et al. [3J, we assume 
a transition from B-DNA to ssDNA. Hence the part of the curve after the transition is fitted by assuming that only 
a single strand remains attached to the cantilever, the second free strand being splitted off. One thus finds a good 
fit (within experimental error bars) by keeping the same parameter values as for Fig. [I] but with a smaller k ss = 0.75 
(instead of 1.5 in Fig. [IJ. This difference might be due to the different base-pair sequence. Two important remarks 
can be done. First, the ratio 7 = a ss /aB is around 2.85, which is significantly larger than the geometrical expected 
value of 2.1. It can be attributable either to the use of Eq. (12 1 with a small Kuhn length, or to the fact that the 
z = reference was not well set in the experiment. Another explanation could be that the spontaneous curvature of 
the AT sequence [M] decreases the effective bp length in the B state. 

Second and more importantly, we are not able to properly fit the transition, especially the part of the curve which 
is after the transition using Eq. (47) (black solid curve). This is due to the fact that, at the transition, the number of 
degrees of freedom increases since for ssDNA the effective bond length is a = 0.285 a ss . Taking into account this fact 
would require drastically a different model. Hence we conclude that the B to ssDNA transition cannot be explained 
by such a mesoscopic model where the monomeric unit remains unchanged through the transition. Writing a model 
where the number of monomers ./V is force-dependent remains challenging. 



CONCLUSION 



This work intends to clarify some issues related to the mesoscopic modeling of DNA molecules subject to an external 
force, which can be applied, for example, by an optical tweezer or an AFM tip. We have focused on both ssDNA 
and dsDNA molecules because it has been suggested that dsDNA can denaturate under load, thus leading to two 
unpaired, single strands. 



We have first addressed the modeling of ssDNA under load and provide an analytical formula, Eq. (16), that allows 



to obtain a very good fit on a wide range of forces, from to 1 nN. Our main conclusion is that this molecule cannot 
be accurately modeled by a polymer, the monomer of which is a nucleobase. This finding is remarkable: whereas 
mesoscopic DNA models using the nucleobase as an elementary building block are usually relevant at small forces, 
the strong force regime requires the use of smaller, sub-nucleobase monomers. Under strong load, the constitutive 
chemical elements of a nucleobase play a significant role because they do not constitute a perfectly rigid entity. This 
is particularly important if one wishes to model the B to ss or S to ss transition: writing a mesoscopic model where 
the nature of the monomers and their number varies continuously with an external parameter (here the applied force) 
is challenging and, to our knowledge, has never been undertaken. 

As far as dsDNA is concerned, we have proposed a new generalization of the extensible, discrete Marko-Siggia 
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Parameters values: Na B = 0.24 fj,m, k b = 147, 7 = 2.85, k 3S = k bss = 0.75, 6 = 0.285, fi = 5.2, J = 1.5, E B = 800. 




formula 34, 35] to a two-state model, which describes successfully force-extension transitions of semiflexible polymers. 
During their first transition, A-phage or poly(dG-dC) DNAs remain in a duplex form, that we have called the "S 
form" as others, where the helix is unwound and successive bases are unstacked. Eq. (47 1, enables us to fit accurately 
experimental B to S transitions and its validity is corroborated by an exact transfer matrix approach (which is 
computationally more complex). These fits, done on different data sets of A-DNA and poly(dG-dC) (Figs. [2] and [3|, 
yield consistently the same parameter values for the S-DNA bending rigidity Kg ~ 4 and monomer length 0,3/0.3 — 
1.7 — 1.9, and the Ising parameters, fj, = 4.5 and J = 1.7 for poly(dG-dC), and /1 ~ 4, J = 2 for A-DNA. The latter 
are consistent with the small but non-negligible influence of the base sequence. 

In contrast, this model is not able to fit the transition observed for poly(dA-dT) (see Fig. [5]). Our conclusion is that 
poly(dA-dT) is subject to a peeling transition where dsDNA is denaturated, thus confirming previous analysis [2]. 
In this ss form, a single strand seems to remain under load. While the B to S transition can be accurately modeled 
because the basic unit (a base) of the model remains the same in the B and S forms, this is not the case for a B to 
ss transition. Note that for the S to ss transition observed for poly(dG-dC) at 320 pN, entropic effects are negligible 
and our model still yields an acceptable fit (see Fig. |4]Jb), using the same parameters for ssDNA as in Fig. [TJ but not 
for A-DNA. It must be emphasized that to conclude these points it was central to be able to fit ssDNA stretching 
curves on the nanoNewton range. 

Put together, these results suggest that when increasing the strength of the base-pairing between both strands (at 
a given salt concentration), the nature of the transition changes. At low base-pairing strength, e.g. for poly(dA-dT), 
unpairing occurs at sufficient low forces so that both unpairing and unstacking transitions are simultaneous |14j . At 
high enough base-pairing strength, in A-phage or poly(dA-dT), unstacking occurs first, leading to the S form, and it 
is followed at stronger forces by unpairing, leading to the ss form. The effect of base sequence is much smaller for the 
B to S transition, maybe because it is essentially the base stacking which is modified during the transition. 

In a very recent work, Zhang et al. |19] observed, using magnetic tweezers, that B-DNA has two different structural 
transitions around 60-70 pN, selected by the temperature or the salt concentration. They differ thermodynamically by 
the sign of df c /dT at the transition, slightly positive for the nonhysteretic (probably B to S) transition and positive 
for the hysteretic (B to ss) one. Our approach yields a negative slope for the B to S transition in the accessible 
temperature window, as seen in Fig. [4ja) , which seems contradictory. However, in our mesoscopic model, we did 
not consider solvent or counterions entropy which are implicitly included in our parameter fj,. Taken /i as a function 
of temperature might allow us to reconcile these two approaches. This work is in progress. Our result, Eq. (47 1, 



may thus be useful to study, in systematic experiments, the role of the variation of salt concentrations [52 and base- 
sequence on stretching transitions. Since our study aims to bridge the gap between force-extension curves and thermal 
denaturation profiles, one thus could benefit from the huge quantity of works done in DNA melting [49 ] 150 ] 153] . 

Finally, we use equilibrium models for describing the transitions and do not consider effects of loading rates [2J [50] 
[S3] , and hysteresis [T7HIII] • Note however that it has been shown that rehybridisation [55] of one strand or the closure 
of one denaturation bubble are very long processes of several /is depending on the DNA length, and one can expect 
that such equilibrium approaches remain valid at moderate loading rates. 
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Appendix A: How to fit ssDNA experimental force-extension curves 

We recall the equation derived in the text for the force-extension curve of ssDNA, f(z): 



fc B T L y niyj " \ 1 + u(k) V1 + 4k 2 ; Y [1 - «(1 + U a i(f))/L 

where the ssDNA contour length L, the bending rigidity modulus n (in k-gT units), and the "effective" monomer size a 
(which can differ from the nucleobase length a ss ~ 0.7 nm) are the three unknown parameters. The Langevin function 
is u(x) = cothz - 1/x and the non-linear stretching polynomial is U n i(f) = 1.172777 / - 3.731836 / 2 + 4.118249 f 3 
where / must be in units of 10 nN. 

Since Eq. ( Al | is highly non-linear, a convenient and extremely accurate simplification is to use / defined by: 



= -( V- u s*; - . 1 ) + J , \ v/rr^ (A2) 

k B T L\1 + u(k) V1 + 4«V V^-z/L) 2 



in the argument of U n \ in Eq. ( Al ). Computing f(z) is then very simple using, for instance, a spreadsheet. 



Appendix B: How to fit B to S-DNA transition in experimental force-extension curves 



The fitting procedure is the following: (i) We fit the low force regime, < / < 20 pN, using Eq. |T3| ) (identical 
to Eq. ( A2 ) above) by fixing the known B-DNA bending modulus kb = 147 [31] (32J • This amounts to fitting the 
scale-parameter for the z axis, the contour length of B-DNA Nas (this step is done here only for Fig. |2|a) since 
the data have already been rescaled in the z-scale for the other data sets), (ii) The stretching modulus Eb is then 
determined by fitting the whole data before the transition (0 < / < 60 pN). (hi) The ratio 7 = as /as and the 
bending modulus of the S form, Ks, are fixed by fitting only the data after the transition using Eq. ( A2 1. (iv) Finally 
the transition is fitted using Eq. (47) which is: 
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(B2) 
(B3) 
(B4) 

(B5) 
(B6) 



Since Na B , 7, k b , ks and E B are known thanks to steps (i)-(iii), this last step yields the values of the parameters 
jj, and J, thus fixing the position (defined by (Xq = 0) and the width of the transition respectively. We have checked 
that choosing ks = k B s does not change significantly the results. 
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